#---------------------
#Release Date: 20191208
#---------------------

#---------------------------------------------------------
proc Mode_Proc { ModeNum } {
	set lambda [eigen -fullGenLapack $ModeNum]
	Mode_Output_Proc  $lambda "Period"
}

#---------------------------------------------------------
proc Mode_Output_Proc { lambda fileName } {
	set pi 3.1415926
	set period "$fileName.txt"
	set Periods [open $period "w"]
	set i 1
	foreach lam $lambda {
		set period [expr 2*$pi / sqrt($lam)]
		set period [format "%.3f" $period]
		set str "Mode "; append str $i; append str ": "; append str $period
		puts $Periods "$str"; puts "$str"
		incr i
	}
	close $Periods
	record
}

#---------------------------------------------------------
proc Gravity_Proc { Step } {
	set incr [expr 1./$Step]	
	constraints Transformation
	numberer RCM
	system UmfPack
	test EnergyIncr 1.0e-6 200
	integrator LoadControl $incr
	algorithm Newton
	analysis Static
	analyze $Step
	puts "Gravity Done."
	loadConst -time 0.0
}

#---------------------------------------------------------
proc Rayleigh_Proc { Mode1 Mode2 Damp } {
	set lambdaN [eigen -fullGenLapack [expr $Mode2]]
	Mode_Output_Proc $lambdaN "Rayleigh"
	set lambdaI [lindex $lambdaN [expr $Mode1-1]]
	set lambdaJ [lindex $lambdaN [expr $Mode2-1]]
	set omegaI [expr pow($lambdaI,0.5)]
	set omegaJ [expr pow($lambdaJ,0.5)]
	set alphaM [expr $Damp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]
	set betaKcurr [expr 2.*$Damp/($omegaI+$omegaJ)]
	rayleigh $alphaM $betaKcurr 0 0
	puts "Rayleigh Damping Defined"
}

#---------------------------------------------------------
proc MRayleigh_Proc { Mode1 Mode2 Damp } {
	set pi 3.1415926
	set omegaI [expr 2*$pi / $Mode1]
	set omegaJ [expr 2*$pi / $Mode2]
	set alphaM [expr $Damp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]
	set betaKcurr [expr 2.*$Damp/($omegaI+$omegaJ)]
	rayleigh $alphaM $betaKcurr 0 0
	puts "Rayleigh Damping Done: Period_1 $Mode1, Period_2 $Mode2."
}

#---------------------------------------------------------
proc GMDefine_Proc { LoadTag GMFileList GMDirList GMFactList dt } {
	foreach GMFile $GMFileList GMDir $GMDirList GMFact $GMFactList {
		#File Exist Then Continue (Ensure All)
		if { [file exists $GMFile] == 0} {
			puts "Error: $GMFile doesn't Existed!"
			return -1
		}
		#Define Series
		set AccelSeries "Series -dt $dt -filePath $GMFile -factor $GMFact"
		#Define Load Case
		pattern UniformExcitation $LoadTag $GMDir -accel "Series -dt $dt -filePath $GMFile -factor $GMFact"
		incr LoadTag
	}
	return 0
}

#---------------------------------------------------------
proc GMDAnalyze_Proc { Duration dt { Name "" } } {
	constraints Transformation
	numberer RCM
	system UmfPack
	test NormDispIncr 0.05 10
	algorithm NewtonLineSearch -tol 0.1
	integrator Newmark 0.55 0.2765625 
	analysis Transient
	set Steps [expr int($Duration / $dt)]
	set Process 5
	puts "*$Name 0%"
	for { set step 1 } { $step <= $Steps } { incr step } {
		set ok [ GMDVariableAnalyze $dt ]
		if {$ok != 0} {
			return -1
		}
		set accu_process [expr $step * 100 / $Steps ]
		if { $accu_process >= $Process} {
			puts "*$Name $Process%"
			set Process [expr $Process + 5]
		}
	}
	return $Duration
}


#---------------------------------------------------------
proc GMDVariableAnalyze { dtMax { numSublevels 10 } { iterMin 10 } } {
	#Current Level
	set CurLevel 1
	#Accumulate time
	set AccTime 0
	#Unitl Time Beyond dt
	while { $AccTime < $dtMax } {
		if { $CurLevel <= 5 } {
			test NormDispIncr 0.05 [expr $iterMin * $CurLevel]
		} else {
			test RelativeNormDispIncr 0.05 [expr $iterMin * $CurLevel]
		}
		#Analysis
		set ok [analyze 1 [expr $dtMax / $CurLevel]]
		#Success or not
		if { $ok == 0 } {
			#Get Accumulate Analysis Time
			set AccTime [expr $AccTime + [expr $dtMax / $CurLevel ]]
		} else {
			#SubLevel Analysis
			incr CurLevel
			puts "Try SubLevel: $CurLevel"
			#Check Level OutRange
			if { $CurLevel > $numSublevels } {
				return -1
			}
		}
	}
	test NormDispIncr 0.05 10
	return 0
}
#---------------------------------------------------------
proc VibrationAnalysis { Duration dt nodeTag dof } {
	#using Math.h
	namespace import ::tcl::mathfunc::*
	#MaxAbsDisp
	set max_disp 0
	#Steps Number
	set Steps [expr int($Duration / $dt)]
	#Analysis Step by Step
	for { set step 1 } { $step <= $Steps } { incr step } {
		#Vibration Analysis
		set ok [ GMDVariableAnalyze $dt ]
		#Fail or Not
		if {$ok != 0} {
			return -999
		}
		#Get MaxDisplacement
		set max_disp [get_Max $max_disp [abs [nodeDisp $nodeTag $dof]]]
	}
	return $max_disp
}

#---------------------------------------------------------
proc VibrationControlAnalysis { Duration dt nodeTag dof tol { Name "" } {MaxTimes 10} } {
	#using Math.h
	namespace import ::tcl::mathfunc::*
	#Duration Must Bigger Than 1
	if { $Duration <= 1. } {
		set Duration 1.
	}
	#Strat VibationAnalysis
	puts "*$Name SelfVibration Analysis..."
	#Vibration Times
	set times 2
	#i MaxDisplacement
	set MaxDispi [VibrationAnalysis $Duration $dt $nodeTag $dof]
	#j MaxDisp
	set MaxDispj [VibrationAnalysis $Duration $dt $nodeTag $dof]
	#Check AnalysisState
	if { $MaxDispj == -999 } {
		return -1
	}
	#Relative Error
	set Error [expr [abs [expr $MaxDispi - $MaxDispj]] / $MaxDispi]
	#While The Displacement Satisfy the Tolerance
	while { $Error >= $tol  } {
		#MaxDisplacement
		set MaxDispi [expr $MaxDispj]
		set MaxDispj [VibrationAnalysis $Duration $dt $nodeTag $dof]
		#Check AnalysisState
		if { $MaxDispj == -999 } {
			return -1
		}
		#Relative Error
		set Error [expr [abs [expr $MaxDispi - $MaxDispj]] / $MaxDispi]
		#VibationAnalysis Times
		incr times
		#Check Time Beyond Analysis Limit Time
		if { $times >= $MaxTimes } {
			break
		}
	}
	#Addictional Analysis Times
	set AdditTime [expr $Duration * $times]
	return $AdditTime
}

#---------------------------------------------------------
proc GetFileMaxAbsValue { GMFile } {
	#using Math.h
	namespace import ::tcl::mathfunc::*
	#File Exist Then Continue (Ensure All)
	if { [file exists $GMFile] == 0} {
			puts "Error: $GMFile doesn't Existed!"
			return -1
	} else {
		#initial value
		set max_abs_value 0
		#Read File
		set fp [open $GMFile r]
		#Read line by line
		while { [gets $fp data] >= 0 } {
			set value [abs $data]
			set max_abs_value [get_Max $value $max_abs_value]
		}
		#Close File
		close $fp
		#RETURN
		return $max_abs_value
	}
}

## Function to Get the Bigger Value
proc get_Max {value_1 value_2} {
    if { $value_1 > $value_2 } {
        return $value_1
    } else {
        return $value_2
    }
}

#---------------------------------------------------------
proc GetFilePrtNum { GMFile } {
	#File Exist Then Continue (Ensure All)
	if { [file exists $GMFile] == 0} {
			puts "Error: $GMFile doesn't Existed!"
			return -1
	} else {
		#initial value
		set prtnum 0
		#Read File
		set fp [open $GMFile r]
		#Read line by line
		while { [gets $fp data] >= 0 } {
			incr prtnum
		}
		#Close File
		close $fp
		#RETURN
		return $prtnum
	}
}










